function matrix_utils()
% MATRIX_UTILS 矩阵计算工具函数集合
% 
% 提供各种矩阵计算中常用的工具函数

fprintf('=== 矩阵计算工具函数 ===\n\n');
fprintf('本文件包含以下工具函数:\n');
fprintf('1. generate_test_matrices - 生成测试矩阵\n');
fprintf('2. matrix_properties - 分析矩阵性质\n');
fprintf('3. condition_analysis - 条件数分析\n');
fprintf('4. timing_utils - 计时工具\n');
fprintf('5. error_analysis - 误差分析\n');
fprintf('6. visualization_utils - 可视化工具\n\n');

end

function A = generate_random_matrix(m, n, condition_number)
% GENERATE_RANDOM_MATRIX 生成指定条件数的随机矩阵
% 
% 输入:
%   m, n - 矩阵维数
%   condition_number - 期望的条件数
% 输出:
%   A - m×n 随机矩阵

if nargin < 3
    condition_number = 1;
end

% 生成随机正交矩阵
[U, ~] = qr(randn(m, m));
[V, ~] = qr(randn(n, n));

% 生成奇异值
sigma_max = 1;
sigma_min = sigma_max / condition_number;
sigma = logspace(log10(sigma_min), log10(sigma_max), min(m, n));

% 构造矩阵
S = zeros(m, n);
for i = 1:min(m, n)
    S(i, i) = sigma(i);
end

A = U * S * V';

end

function A = generate_symmetric_matrix(n, condition_number, definite_type)
% GENERATE_SYMMETRIC_MATRIX 生成对称矩阵
% 
% 输入:
%   n - 矩阵维数
%   condition_number - 条件数
%   definite_type - 'positive', 'negative', 'indefinite'

if nargin < 3
    definite_type = 'positive';
end

% 生成随机正交矩阵
[Q, ~] = qr(randn(n, n));

% 生成特征值
switch definite_type
    case 'positive'
        lambda_max = 1;
        lambda_min = lambda_max / condition_number;
        eigenvals = logspace(log10(lambda_min), log10(lambda_max), n);
    case 'negative'
        lambda_min = -1;
        lambda_max = lambda_min / condition_number;
        eigenvals = -logspace(log10(-lambda_max), log10(-lambda_min), n);
    case 'indefinite'
        pos_vals = logspace(log10(1/condition_number), 0, floor(n/2));
        neg_vals = -logspace(log10(1/condition_number), 0, ceil(n/2));
        eigenvals = [pos_vals, neg_vals];
        eigenvals = eigenvals(randperm(n));
end

A = Q * diag(eigenvals) * Q';

end

function A = generate_banded_matrix(n, bandwidth, condition_number)
% GENERATE_BANDED_MATRIX 生成带状矩阵
% 
% 输入:
%   n - 矩阵维数
%   bandwidth - 带宽
%   condition_number - 条件数

A = zeros(n, n);

% 生成随机带状矩阵
for i = 1:n
    for j = max(1, i-bandwidth):min(n, i+bandwidth)
        A(i, j) = randn();
    end
end

% 调整条件数
if nargin >= 3 && condition_number > 1
    % 通过SVD调整条件数
    [U, S, V] = svd(A);
    s = diag(S);
    s_new = logspace(log10(max(s)/condition_number), log10(max(s)), length(s));
    A = U * diag(s_new) * V';
    
    % 保持带状结构
    for i = 1:n
        for j = 1:n
            if abs(i-j) > bandwidth
                A(i, j) = 0;
            end
        end
    end
end

end

function props = analyze_matrix_properties(A)
% ANALYZE_MATRIX_PROPERTIES 分析矩阵性质
% 
% 输入: A - 矩阵
% 输出: props - 包含矩阵性质的结构体

props = struct();

[m, n] = size(A);
props.size = [m, n];
props.is_square = (m == n);
props.rank = rank(A);
props.is_full_rank = (props.rank == min(m, n));

if props.is_square
    props.determinant = det(A);
    props.trace = trace(A);
    props.is_singular = (abs(props.determinant) < eps);
    props.condition_number = cond(A);
    
    % 对称性检查
    props.is_symmetric = norm(A - A', 'fro') < eps * norm(A, 'fro');
    
    if props.is_symmetric
        eigenvals = eig(A);
        props.eigenvalues = eigenvals;
        props.is_positive_definite = all(eigenvals > eps);
        props.is_negative_definite = all(eigenvals < -eps);
        props.is_positive_semidefinite = all(eigenvals >= -eps);
        props.is_indefinite = any(eigenvals > eps) && any(eigenvals < -eps);
    end
    
    % 特殊矩阵类型检查
    props.is_diagonal = norm(A - diag(diag(A)), 'fro') < eps * norm(A, 'fro');
    props.is_upper_triangular = norm(tril(A, -1), 'fro') < eps * norm(A, 'fro');
    props.is_lower_triangular = norm(triu(A, 1), 'fro') < eps * norm(A, 'fro');
    
    % 正交性检查
    if ~props.is_singular
        props.is_orthogonal = norm(A * A' - eye(n), 'fro') < eps * n;
    else
        props.is_orthogonal = false;
    end
end

% 范数计算
props.norm_1 = norm(A, 1);
props.norm_2 = norm(A, 2);
props.norm_inf = norm(A, inf);
props.norm_fro = norm(A, 'fro');

end

function display_matrix_properties(A, name)
% DISPLAY_MATRIX_PROPERTIES 显示矩阵性质
% 
% 输入: A - 矩阵, name - 矩阵名称

if nargin < 2
    name = 'A';
end

props = analyze_matrix_properties(A);

fprintf('=== 矩阵 %s 的性质 ===\n', name);
fprintf('维数: %d×%d\n', props.size);
fprintf('秩: %d\n', props.rank);
fprintf('满秩: %s\n', mat2str(props.is_full_rank));

if props.is_square
    fprintf('行列式: %.6e\n', props.determinant);
    fprintf('迹: %.6f\n', props.trace);
    fprintf('条件数: %.2e\n', props.condition_number);
    fprintf('对称: %s\n', mat2str(props.is_symmetric));
    
    if props.is_symmetric
        fprintf('正定: %s\n', mat2str(props.is_positive_definite));
        fprintf('半正定: %s\n', mat2str(props.is_positive_semidefinite));
        fprintf('不定: %s\n', mat2str(props.is_indefinite));
    end
    
    fprintf('对角: %s\n', mat2str(props.is_diagonal));
    fprintf('上三角: %s\n', mat2str(props.is_upper_triangular));
    fprintf('下三角: %s\n', mat2str(props.is_lower_triangular));
    fprintf('正交: %s\n', mat2str(props.is_orthogonal));
end

fprintf('范数:\n');
fprintf('  ||A||₁ = %.6f\n', props.norm_1);
fprintf('  ||A||₂ = %.6f\n', props.norm_2);
fprintf('  ||A||∞ = %.6f\n', props.norm_inf);
fprintf('  ||A||F = %.6f\n', props.norm_fro);
fprintf('\n');

end

function results = benchmark_algorithms(algorithms, test_matrices, varargin)
% BENCHMARK_ALGORITHMS 算法性能基准测试
% 
% 输入:
%   algorithms - 算法函数句柄的cell数组
%   test_matrices - 测试矩阵的cell数组
%   varargin - 可选参数

p = inputParser;
addParameter(p, 'num_runs', 5, @isnumeric);
addParameter(p, 'display_results', true, @islogical);
parse(p, varargin{:});

num_runs = p.Results.num_runs;
display_results = p.Results.display_results;

num_algs = length(algorithms);
num_matrices = length(test_matrices);

results = struct();
results.times = zeros(num_algs, num_matrices);
results.errors = zeros(num_algs, num_matrices);
results.algorithm_names = cell(num_algs, 1);
results.matrix_names = cell(num_matrices, 1);

% 获取算法名称
for i = 1:num_algs
    func_info = functions(algorithms{i});
    results.algorithm_names{i} = func_info.function;
end

% 运行基准测试
for j = 1:num_matrices
    A = test_matrices{j};
    results.matrix_names{j} = sprintf('Matrix_%d (%d×%d)', j, size(A));
    
    for i = 1:num_algs
        times = zeros(num_runs, 1);
        
        for run = 1:num_runs
            tic;
            try
                result = algorithms{i}(A);
                times(run) = toc;
                
                % 计算误差（如果可能）
                if isstruct(result) && isfield(result, 'error')
                    results.errors(i, j) = result.error;
                end
            catch ME
                times(run) = inf;
                results.errors(i, j) = inf;
                if display_results
                    fprintf('算法 %s 在矩阵 %d 上失败: %s\n', ...
                            results.algorithm_names{i}, j, ME.message);
                end
            end
        end
        
        results.times(i, j) = median(times);
    end
end

% 显示结果
if display_results
    fprintf('\n=== 算法性能基准测试结果 ===\n');
    fprintf('时间 (秒):\n');
    fprintf('%15s', 'Algorithm');
    for j = 1:num_matrices
        fprintf('%12s', sprintf('Mat_%d', j));
    end
    fprintf('\n');
    
    for i = 1:num_algs
        fprintf('%15s', results.algorithm_names{i});
        for j = 1:num_matrices
            if isfinite(results.times(i, j))
                fprintf('%12.4f', results.times(i, j));
            else
                fprintf('%12s', 'FAIL');
            end
        end
        fprintf('\n');
    end
end

end

function plot_convergence(errors, varargin)
% PLOT_CONVERGENCE 绘制收敛曲线
% 
% 输入:
%   errors - 误差向量或矩阵
%   varargin - 可选参数

p = inputParser;
addParameter(p, 'semilogy', true, @islogical);
addParameter(p, 'title', '收敛曲线', @ischar);
addParameter(p, 'xlabel', '迭代次数', @ischar);
addParameter(p, 'ylabel', '误差', @ischar);
addParameter(p, 'legend', {}, @iscell);
parse(p, varargin{:});

figure('Name', p.Results.title);

if p.Results.semilogy
    semilogy(errors, 'LineWidth', 2);
else
    plot(errors, 'LineWidth', 2);
end

xlabel(p.Results.xlabel);
ylabel(p.Results.ylabel);
title(p.Results.title);
grid on;

if ~isempty(p.Results.legend)
    legend(p.Results.legend, 'Location', 'best');
end

end

function plot_matrix_structure(A, varargin)
% PLOT_MATRIX_STRUCTURE 可视化矩阵结构
% 
% 输入: A - 矩阵

p = inputParser;
addParameter(p, 'title', '矩阵结构', @ischar);
addParameter(p, 'colormap', 'default', @ischar);
parse(p, varargin{:});

figure('Name', p.Results.title);

% 稀疏矩阵用spy，稠密矩阵用imagesc
if issparse(A) || nnz(A) / numel(A) < 0.1
    spy(A);
    title([p.Results.title, ' (稀疏模式)']);
else
    imagesc(abs(A));
    colorbar;
    title([p.Results.title, ' (绝对值)']);
    
    if ~strcmp(p.Results.colormap, 'default')
        colormap(p.Results.colormap);
    end
end

axis equal tight;

end

function save_results(results, filename)
% SAVE_RESULTS 保存计算结果
% 
% 输入:
%   results - 结果结构体
%   filename - 文件名

if nargin < 2
    filename = sprintf('matrix_results_%s.mat', datestr(now, 'yyyymmdd_HHMMSS'));
end

save(filename, 'results');
fprintf('结果已保存到文件: %s\n', filename);

end

function results = load_results(filename)
% LOAD_RESULTS 加载计算结果
% 
% 输入: filename - 文件名
% 输出: results - 结果结构体

if exist(filename, 'file')
    data = load(filename);
    results = data.results;
    fprintf('结果已从文件加载: %s\n', filename);
else
    error('文件不存在: %s', filename);
end

end